ONE-seq: epitranscriptome and gene-specific profiling of NAD-capped RNA

Abstract The hub metabolite, nicotinamide adenine dinucleotide (NAD), can be used as an initiating nucleotide in RNA synthesis to result in NAD-capped RNAs (NAD-RNA). Since NAD has been heightened as one of the most essential modulators in aging and various age-related diseases, its attachment to RNA might indicate a yet-to-be discovered mechanism that impacts adult life-course. However, the unknown identity of NAD-linked RNAs in adult and aging tissues has hindered functional studies. Here, we introduce ONE-seq method to identify the RNA transcripts that contain NAD cap. ONE-seq has been optimized to use only one-step chemo-enzymatic biotinylation, followed by streptavidin capture and the nudix phosphohydrolase NudC-catalyzed elution, to specifically recover NAD-capped RNAs for epitranscriptome and gene-specific analyses. Using ONE-seq, we discover more than a thousand of previously unknown NAD-RNAs in the mouse liver and reveal epitranscriptome-wide dynamics of NAD-RNAs with age. ONE-seq empowers the identification of NAD-capped RNAs that are responsive to distinct physiological states, facilitating functional investigation into this modification.


INTRODUCTION
In eukaryotes, 5 ,5 -triphosphate-linked 7-methylguanosine (m 7 G) is the predominant 5 -end cap structure of RNA (m 7 Gppp-RNA or m 7 G-RNA), essential for RNA stability, polyadenylation, splicing, localization, and translation (1)(2)(3). Recently, NAD, the adenine nucleotide containing metabolite, emerged as a non-canonical initiating nucleotide (NCIN) incorporating at the 5 -terminus of RNA to result in NAD-capped RNAs (NAD-RNA), which are estimated to make up more than 0.6% and 1.3% of the mRNA transcripts expressed in the entire transcriptome from mouse liver and kidney, respectively (4,5). NAD capping may define a yet-to-be understood epitranscriptomic mechanism.
Mounting evidence associates aging with transcriptional and metabolic alterations (6)(7)(8), yet how these processes are integrated into the regulation of aging is only beginning to be revealed. NAD is the hub metabolite and redox regent for cells, involving in a wide range of biological processes (9). In rodents and humans, studies have revealed that the NAD level declines with age in critical organs including liver (10). Boosting NAD biosynthesis, on the other hand, extends lifespan in yeast (11), worms (12), flies (13) and mice (14). Given the dynamics of NAD and gene expression over the course of adult lifespan, it is important to better understand how NAD-RNAs are modulated with age and its consequent impact on the progression of aging. Knowing the identity of NAD-modified RNAs is prerequisite for functional studies. However, investigating functional insight of NAD-RNAs has been largely limited by the analytical methods available.
The currently reported NAD-RNA identification methods involve the use of multiple chemo-enzymatic reactions (15). NAD captureSeq was initially reported, utilizing an enzymatic reaction by adenosine diphosphate ribosyl-Animal experimentation C57BL/6 mice were housed in an environmentally controlled room under a 12:12 h light/dark cycle at 23 • C and were fed with commercial mouse chow food and water ad libitum. To examine the efficacy of NudC-catalyzed NAD-RNA elution, livers were dissected from 12-month old C57BL/6 mice. To profile NAD-RNA epitranscriptomes, livers were dissected from 2-and 18-month old C57BL/6 mice. Dissected livers were immediately frozen in liquid nitrogen.

HPLC and LC-MS analysis
Samples were analyzed by an Agilent Technologies 1260 Infinity with a C18 column (Phenomenex Kinetex, 5 m EVO C18, 100Å, 150 × 4.6 mm) monitoring at 260 nm. The mobile phase was composed of (A) 100 mM NH 4 HCO 3 in water and (B) acetonitrile at a flow rate of 0.4 ml/min. The injection volume was 2 l. The gradient was as follows: 0 min 0% B, 2 min 0% B, 4 min 20% B, 17 min 40% B, 19 min 100% B, 24 min 100% B, and 30 min 0% B. The retention time of NAD and ADPR-Biotin (product) was 9.393 min and 9.929 min, respectively. The corresponding compounds were purified from HPLC and then analyzed by LC-MS by Agilent Technologies 6120 Quadrupole with a C18 column (Agilent Poroshell, 120 EC-C18, 2.7 m, 3.0 × 50 mm). The mobile phase was composed of (A) 0.1% formic acid in water and (B) acetonitrile at a flow rate of 0.4 ml/min. The injection volume was 5 l. The gradient was as follows: 0 min 10% B, 4 min 100% B, and 7 min 100% B.

HEEB reaction with NAD-RNA (38 nt)
NAD-RNA (1 g) was reacted with 40 mM and 100 mM HEEB (1 M stock in DMSO) under the catalysis of AD-PRC (25 g/ml) with 0.4 U/l of RNase Inhibitor (Takara Bio, catalog: 2313A) in 100 l of ADPRC reaction buffer (50 mM Na-HEPES pH 7.0, 5 mM MgCl 2 ) at 37 • C for 1 h. Product was purified with Zymo column (RNA Clean & Concentrator-5, Zymo Research, catalog: R1016) according to the instruction of manufacturer and analyzed by an 8% polyacrylamide TBE urea gel. Gel was stained by SYBR Gold (Invitrogen, catalog: S11494) and fluorescence was detected by Typhoon FLA 7000 fluorescent image analyzer (GE Life Science).

Enrichment of biotinylated NAD-RNA (38 nt) and m 7 Gppp-RNA (38 nt)
Magnetic streptavidin beads (6 l, MedChemExpress, catalog: HY-K0208) were washed 3 times with immobilization buffer (10 mM Na-HEPES, 1 M NaCl, 5 mM EDTA) (17), and then incubated with the purified RNA product from HEEB reaction in immobilization buffer, together with 0.4 U/l of RNase Inhibitor, in a thermomixer at room temperature for 30 min. After the beads were washed 5 times with streptavidin wash buffer (50 mM Tris-HCl, pH 7.4, 8 M urea) (17), the biotinylated RNA on the beads was extracted with Trizol LS (Ambion, catalog: 10296010) and chloroform. The upper aqueous layer was further purified with Zymo column (Zymo Research, catalog: R1016). The purified RNA product was analyzed by an 8% polyacrylamide TBE urea PAGE gel. Gel was stained by SYBR Gold and fluorescence was detected by Typhoon FLA 7000 fluorescent image analyzer (GE Life Science).

Dot blot analysis
500 ng RNA was blotted on Hybond N + membrane (GE Healthcare, catalog: RPN203B) which was then crosslinked by UV 254nm (0.18 J/cm 2 ) twice. Membrane was blocked in 5% BSA in PBST (0.1% Tween20 in PBS) for 1 h, followed by incubation with Alexa Fluor® 790 Streptavidin (Jackson, catalog: 016-650-084, 1:10 000 in PBST) at room temperature for 2 h in the dark. Membrane was washed with PBST and PBS for 3 times, respectively, and imaged on the Odyssey LiCor CLx scanner (Li-Cor Biosciences) with the software set to auto-detect the signal intensity at 800 nm channel. Finally, the membrane was stained with methylene blue solution (0.3% w/v methylene blue + 30% v/v ethanol + 70% v/v H 2 O) (19) to visualize the presence of RNA. buffer (100 mM NaCl, 50 mM Tris-HCl pH 7.9, 10 mM MgCl 2 , 100 g/ml Recombinant Albumin) at 37 • C for 30 min. Product was purified with Trizol LS (Ambion, catalog: 10296010) and analyzed by an 8% polyacrylamide TBEurea gel. Gel was stained by SYBR Gold (Invitrogen, catalog: S11494) and fluorescence was detected by Typhoon FLA 7000 fluorescent image analyzer (GE Life Science).

HEEB reaction with total RNA extract
Total RNA was prepared from mouse liver tissues, in accordance with manufacturer's instruction (Takara Bio, catalog: 9108). Total RNAs (100 g) was incubated with 100 mM HEEB (1 M stock in DMSO), ADPRC (25 g/ml, Sigma-Aldrich, catalog: A9106) and 0.4 U/l of RNase Inhibitor (Takara Bio, catalog: 2313B) in 100 l of ADPRC reaction buffer at 37 • C for 1 h. 100 l of DEPC-treated H 2 O was then added and acid phenol/ether extraction was performed to stop the reaction (17). RNAs were precipitated by ethanol, re-dissolved in 100 l of DEPC-treated H 2 O. 5 l of biotinylated RNAs were kept as input.

NudC-catalyzed NAD-RNA elution
After HEEB reaction, biotinylated RNAs were incubated with streptavidin bead particles (6 l, MedChemExpress, catalog: HY-K0208) and 0.4 U/l of RNase Inhibitor (Takara Bio, catalog: 2313B) at 25 • C for 30 min. Beads were washed four times with streptavidin wash buffer (50 mM Tris-HCl (pH 7.4) and 8 M urea), and three times with DEPC-treated H 2 O. To ensure complete elution, biotinconjugated RNAs were replaced from streptavidin beads by incubating with 1 mM biotin buffer (20 l, Sigma-Aldrich, catalog: B4639) at 94 • C for 8 min, followed by incubation with 500 nM NudC (New England Biolabs, catalog: M0607S) in 25 l of NudC reaction buffer (100 mM NaCl, 50 mM Tris-HCl pH 7.9, 10 mM MgCl 2 , 100 g/ml Recombinant Albumin) at 37 • C for 30 min. After NudC treatment, biotinylated-RNAs that are resistant to NudC catalysis, potentially derived from contaminating m 7 G-RNAs, were retained on beads by incubation with high-capacity streptavidin particle (20 l, Thermo Fisher Scientific, catalog: 20357) at 25 • C for 30 min. Eluted RNAs in the supernatant were used for next step.

PolyA-selected RNA sequencing
Input (see above) and NudC-eluted RNAs from four biological replicates of livers from 2-and 18-month old C57BL/6 mice were used for NGS library construction, in accordance with manufacturer's instructions (mRNA-seq Lib Prep Kit for Illumina, Abclonal, catalog: RK20302). Library quality was assessed by Bioanalyzer 2100 (Agilent, United States), and quantification was performed by qRT-PCR with a reference to a standard library. Libraries were pooled together in equimolar amounts to a final 2 nM concentration and denatured with 0.1 M NaOH (Sigma, catalog: 72068). Libraries were sequenced on the Illumina No-vaSeq 6000 system (paired end; 150 bp).

-End blocking by adaptor ligation
1 g of NAD-RNA (38 nt) or m 7 Gppp-RNA (38 nt) was ligated with 5 M 3 adaptor oligo listed in Supplementary  Table S5, in the presence of 10 U T4 RNA ligase 1 (New England Biolabs, catalog: M0202), 10% of PEG8000, 1 mM ATP and 40U RNaseOUT in 20 l of 1x T4 RNA ligase buffer. Reaction was incubated at 16 • C for 16 h. RNAs were purified by Trizol LS (Ambion, catalog: 10296010) according to the instruction of manufacturer and analyzed by an 8% polyacrylamide TBE urea gel. Gel was stained by SYBR Gold (Invitrogen, catalog: S11494) and fluorescence was detected by Typhoon FLA 7000 fluorescent image analyzer (GE Life Science).

Enrichment of RNA oligonucleotide by boronic acid beads
Boronic acid beads (30 l, Smart-Lifesciences, catalog: SA057005) were washed 3 times with binding buffer (50 mM Na-HEPES, 1 M NaCl, PH 7.8), and then incubated with RNA oligonucleotide in binding buffer in a thermomixer at 37 • C for 2 h. Beads were then washed 3 times with wash buffer (50 mM PBS, pH 7.4, 6 M urea); RNA oligonucleotide bound by the beads was extracted with Trizol LS (Ambion, catalog: 10296010)). Purified RNA product was analyzed by an 8% polyacrylamide TBE-urea PAGE gel. Gel was stained by SYBR Gold and fluorescence was detected by Typhoon FLA 7000 fluorescent image analyzer (GE Life Science).

qRT-PCR analysis of NAD-capped RNAs bound by boronic acid beads
Total RNAs were extracted from three biological replicates of livers from 2-month old C57BL/6 mice as described above. For each replicate, total RNAs (50 g) were mixed with polyadenylated spike-in RNAs (1 ng NAD-RNA, 1 ng m 7 Gppp-RNA, and 1 ng of ppp-RNA). RNAs were incubated with yDcpS (New England Biolabs, catalog: M0463) in 50 l of 1X yDcpS reaction buffer at 37 • C for 1 h, to remove m 7 G cap. After purification with Trizol LS (Ambion, catalog: 10296010), yDcpS-treated RNAs were combined with 100 mol oligo d(T)30VN in 100 l 0.5× SSC buffer (Sangon Biotech, catalog: M0463), and incubated at 80 • C for 5 min followed by 50 • C for 60 min. Then, 1 l RNase H (Thermo Fisher Scientific, catalog: EN0202) and 5.6 l 10× Reaction buffer (Thermo Fisher Scientific, catalog: EN0202) were added. The mixture was incubated at 37 • C for 20 min, and purified with Trizol LS (Ambion, catalog: 10296010). This step aims to promote 3 -end adaptor ligation by removing polyA sequence tract from the endogenous transcripts. Prior to 3 -end ligation, polyA-depleted RNAs were treated with 2 U CIAP (Thermo Fisher Scientific, catalog: 18009019) in the presence of 40 U RNaseOUT at 37 • C for 1 h to remove 5 -terminal phosphate. RNA was extracted with Trizol LS (Ambion, catalog: 10296010). 5 g of CIAP-treated RNA products were ligated with 100 M 3 adaptor oligo listed in Supplementary Table S5, in the presence of 10 U T4 RNA ligase 1 (New England Biolabs, catalog: M0202), 10% of PEG8000, 1 mM ATP and 40 U RNaseOUT in 50 l of 1× T4 RNA ligase buffer. Reaction was incubated at 16 • C for 16 h. RNAs were purified by Trizol LS (Ambion, catalog: 10296010), and re-dissolved in 100 l of binding buffer (50 mM Na-HEPES, 1 M NaCl, PH 7.8). RNAs were incubated with boronic acid beads (30 l, Smart-Lifesciences, catalog: SA057005) at 37 • C for 2 h. Beads were washed three times with wash buffer (50 mM PBS and 6 M urea), and once with DEPC-treated H 2 O. Reverse transcription was performed with adaptorbased primer listed in Supplementary Table S5, by Bey-oRT™ II First-Strand cDNA Synthesis Kit (Beyotime, catalog: D7168M). Meanwhile, 5 g of CIAP-treated RNAs without 3 ligation were incubated with boronic acid beads, reverse transcribed with random hexamer primers, and served as input. After reverse transcription, DNA/RNA mixtures were treated with 1 U RNase H at 37 • C for 20 min. While ppp-RNA (500 nt) was served as the baseline, NAD-RNA (500 nt) was used as a positive control, and m 7 G-RNA (500 nt) was used as a negative control. Primers were listed in Supplementary Table S4. To calculate the fold enrichment of each adaptor + reaction from qPCR data, we normalized the Ct value of the target RNA to the Ct of the ppp-RNA. Fold enrichments were calculated by normalizing adaptor + fraction value ( Ct of target RNA normalized to the ppp-RNA) to nonspecific background (the similar Ct calculation of the input), to yield the Ct value. The linear conversion of this Ct resulted in fold enrichment. Significance was assessed by Student's t-test.

Data analysis
All sequencing reads were processed with Trim Galore (v0.6.6) (20) with the parameters '-nextseq 30 -paired' to remove the adapter sequences (AGATCGGAAGAGC) from NovaSeq-platforms, and reads longer than 20 bp were kept. Reads that passed the quality control procedure were kept and mapped to the Mus musculus genome (GRCm38) using STAR (2.7.6a) (21) with default parameters. Mapped read pairs were counted against Gencode (M23) annotations using featureCounts (v2.0.1) with parameters '-p -B -C' (22). Aligned reads were summarized as exon counts and gene body counts using exon annotation and gene body annotation, respectively. Intron read counts were calculated by subtracting exon counts from gene body counts at genelevel (23). Read alignments were converted to the bigwig files by bamCoverage (24) and were visualized at selected genomic regions by IGV (v2.8.3) (25). Sequencing saturation was assessed by randomly subsampling of the original libraries and analysis of the corresponding changes in gene numbers with more than 10 read counts. The gene body coverage of sequencing reads was calculated using RSeQC (26) by counting the number of reads covering each nucleotide position and all transcripts were scaled to 100 bins.
Genes that had zero count in more than 25% (3 out of 12) sequencing libraries were removed before performing differential analysis. Principle component analysis (PCA)

PAGE 7 OF 17
Nucleic Acids Research, 2023, Vol. 51, No. 2 e12 was performed with R function 'prcomp' using a transformed counts matrix by function 'vst'. To identify NAD-RNA from total RNA-seq data, we used R package DE-Seq2 (v1.30.0) (27) to perform differential analysis. Significance of logarithmic fold changes were determined by a Wald test to approximate P values, and genes passing an independent filtering step were adjusted for multiple testing using the Benjamini-Hochberg procedure to yield a false discovery rate (FDR). NAD-RNAs were defined as fold change of normalized transcript counts ≥2 and FDR <0.05 in NudC-treated samples compared to those in input samples. Gene annotation information, such as chromosome, gene-types, gene-lengths and intron coordinates were retrieved from Gencode (M23) annotations. The violin plot, bar plot, line chart and scatter plot were generated by R package ggplot2 (v3.3.2) (28).
Pathway enrichment analysis was performed using R package gprofiler2 (v0.2.1) (29). Pathways were defined as the molecular pathways of Reactome (BioMart releases: 2021-5-7) (30) and the biological processes (BPs) of GO (BioMart releases: 2021-05-01) (31) without less reliable GO annotations (IEAs) that are not manually reviewed. Size of enriched gene sets were limited to range between 5 and 350. Multiple testing corrections were performed with gprofiler2 built-in 'g SCS' method and terms with adjusted P-value <0.05 were considered as significantly enriched. Top 10 significantly enriched pathways were visualized in bar plot. In addition, the results of the pathway analysis were visualized using the EnrichmentMap App (v1.1.0) in Cytoscape (v3.8.2) (32). Network maps were generated for nodes with FDR <0.05 and nodes sharing gene overlaps with Jaccard coefficient >0.60 were connected by a blue line (edge). Clusters of related pathways were identified and annotated using the AutoAnnotate App (v1.3.5) in Cytoscape that uses a Markov Cluster (MCL) algorithm which connects pathways by shared keyword in the description of each pathway. The resulting clusters of pathways were manually reviewed and were designated as the major pathways in a circle (33).

ONE-seq-based gene specific analysis of NAD capping
Total RNA was extracted from four biological replicates of livers from 2-month old C57BL/6 mice as described above. For each replicate, 100 g total RNA, mixed with 1 ng of NAD-RNA (106 nt) and 1 ng of ppp-RNA (106 nt), was incubated with 100 mM HEEB (1 M stock in DMSO) with ADPRC (25 g/ml) in 100 l of ADPRC reaction buffer (50 mM Na-HEPES pH 7.0, 5 mM MgCl 2 ) at 37 • C for 1 h. 100 l of DEPC-treated H 2 O was then added and acid phenol/ether extraction was performed to stop the reaction (17). RNAs were precipitated by ethanol, re-dissolved in 100 l of DEPC-treated H 2 O. 1 l of biotinylated RNAs were kept as input. After NudC-catalyzed NAD-RNA elution, input and NudC-eluted RNAs were used for reverse transcription. Quantification of transcript abundance by qRT-PCR was performed using SYBR Green master mix (Vazyme, catalog: Q111-02) following the manufacturer's protocol with primers listed in Supplementary Table S4. ppp-RNA (106 nt) was served as the baseline, and NAD-RNA (106 nt) was used as an internal positive control. To calculate the enrichment from qRT-PCR data, the Ct value of the target gene was first normalized to the Ct of the ppp-RNA (baseline). Next, the normalized NudC + fraction value ( Ct of the target gene normalized to the ppp-RNA) was normalized to the background ( Ct calculation for the gene in the input), to yield the Ct value. The linear conversion of this Ct resulted in the fold enrichment.

The workflow of ONE-seq
Compared to previous methods that require multiple reactions, we introduced HEEB (N-[2-(2-hydroxyethoxy) ethyl]biotinamide) as a reactant, allowing only one reaction to biotinylate NAD-capped RNAs. To avoid contaminating signals, we designed a NudC-based post-treatment to elute biotin-conjugated RNAs specifically derived from NAD, but not m 7 G-capped transcripts from streptavidin beads. Therefore, our method circumvented laborious steps to purify mRNAs followed by clearance of m 7 G-capped RNAs.
With significantly simplified procedures, our method enabled NAD-RNA profiling directly from total RNA. Based on the same platform, qRT-PCR analysis on specific NAD-RNAs could be readily applied. We thereby named our method ONE-seq, through Onestep chemo-enzymatic reaction to biotinylate NAD-capped RNA for streptavidin binding, followed by the nudix phosphohydrolase NudC-catalyzed Elution, to specifically harvest NAD-capped RNA from streptavidin beads for highthroughput sequencing (Figure 1). First, HEEB has a terminal hydroxyl group as the nucleophile and a biotin group as the affinity tag. Catalyzed by ADPRC, nicotinamide moiety of NAD can be replaced by HEEB via nucleophilic reaction and simultaneously biotinylated, which subsequently can be enriched by streptavidin beads. Second, NudC, known for its ability to hydrolyze the diphosphate, but not triphosphate, linkage, detaches HEEB-RNAs specifically derived from NAD-capped transcripts from streptavidin beads. At this step, contaminating m 7 G-RNAs that also react with HEEB are retained on the bead. Third, eluted RNAs can be used for epitranscriptome-wide profiling as well as genespecific qRT-PCR analysis.

The feasibility of one-step chemo-enzymatic reaction
We tested the feasibility of one-step chemo-enzymatic reaction. First, we performed an ADPRC catalyzed reaction between NAD molecule and HEEB. HPLC and LC-MS confirmed a product corresponding to the biotinylated NADderived structure (Figure 2A and Supplementary Figure  S1). Second, we subjected a synthetic 38-nucleotide (nt) NAD-capped RNA to the HEEB reaction, resulting in an ADPRC-dependent yielding, evidenced by the accumulation of an upper band in the TBE-UREA gel ( Figure 2B). By incubation with the streptavidin beads, we showed the evidence that only the upper band from the reaction was retained by the streptavidin beads, while the lower band of non-biotinylated form was discarded with flow-through ( Figure 2C). In addition, we applied avidin-conjugated fluorophore to detect biotinylation. Imaging-based dot blot assay corroborated the RNA product being biotin-tagged in ethyl]-biotinamide, in red) has a terminal hydroxyl group as the nucleophile and a biotin group as the affinity tag. In the presence of ADPRC, nicotinamide moiety of NAD can be replaced by HEEB via nucleophilic reaction and simultaneously biotinylated, which subsequently can be enriched by streptavidin beads. NudC (in red) detaches HEEB-RNAs specifically derived from NAD-capped transcripts from streptavidin beads, while contaminating m 7 G-RNAs that also react with HEEB remain retained on the bead. Eluted RNAs can be used for epitranscriptome-wide profiling as well as gene-specific qRT-PCR analysis. the presence of ADPRC ( Figure 2D). To examine the specificity, we subjected a mixture of m 7 Gppp-RNA (45 nt) and NAD-RNA (38 nt) with 100 mM HEEB to the reaction. While NAD-RNA reacted with HEEB, as indicated by the presence of a biotinylated form, m 7 G-RNA had no such a product, reflecting specificity ( Figure 3A). However, we noted that HEEB, only at high concentration (400 mM), reacted with m 7 Gppp-RNA, yielding an upper band in the gel ( Figure 3B). Although we used 100 mM HEEB for the subsequent ONE-seq experiments, we worried that, biotinylated m 7 Gppp-RNA, even at low level, would cause falsepositive signals.

The efficacy of NudC
We evaluated the efficacy of NudC. Our rationale lies at the fact that NudC can cleave the pyrophosphate group of NAD-RNA, while being inactive to the m 7 G cap that comprises triphosphate moiety. As noted, NudC has been previously used by the CapZyme-seq method to study the global profile of NCIN-capped RNAs from total RNA extract (34). However, CapZyme-seq cannot specifically identify NAD-RNAs in vivo (34). Moreover, studies reported non-specific activity of NudC on m 7 G-RNA in the presence of Mn 2+ and upon prolonged incubation (35,36); we therefore performed NudC without Mn 2+ and within 30 min of treatment. We designed different spike-in RNAs to assess the performance of NudC treatment.
First, 38 nt RNA oligos with either NAD or m 7 G-cap were tested. NudC was able to de-cap NAD-RNA (38 nt) but not m 7 Gppp-RNA (38 nt) as shown by a lower-sized band corresponding to the de-capped product ( Figure 3C). We further subjected a mixture of NAD-RNA (38 nt) and m 7 Gppp-RNA (45 nt) to the same reaction. The result clearly showed that NudC was able to selectively de-cap NAD-capped (38 nt) but not m 7 G-capped RNA (45 nt) ( Figure 3D). We tested the ability of NudC to elute biotinylated RNAs bound by streptavidin beads. To ensure complete elution, we first used biotin to replace all bound RNAs from the streptavidin beads, followed by NudC treatment. Then we re-applied high-capacity streptavidin beads to capture biotinylated-RNAs that were resistant to NudC treatment. At this step, contaminating m 7 G-RNAs that remain biotinylated are retained on the beads, while decapped NAD-RNAs are eluted. Treatment by NudC efficiently eluted RNAs derived from NAD-capped (38 nt) but not m 7 Gppp-capped forms (38 nt) from streptavidin beads ( Figure 3E and F). A control experiment corroborated that m 7 Gppp-RNA (38 nt), but not NAD-RNA (38 nt) was sensitive to yDcpS treatment, a decapping enzyme that hydrolyzes the triphosphate linkage of m 7 G-capped RNA (Supplementary Figure S2A).
Third, we synthesized two long RNA spike-ins with identical sequence (500 nt) but had either NAD or m 7 G-cap, followed by polyA tails. Presumably, endogenous transcripts may contain both NAD and m 7 G-capped forms, though the percentage may differ for particular genes. The rationale of this design is to mimic endogenous genes that have either low (0% or 1%) or relatively high (5% or 10%) NAD modification. Total RNAs were mixed with spike-ins that had different ratios of NAD-RNA and m 7 G-RNA. The sample that contained 100% m 7 G-RNA spike-in represented a gene with no NAD capping, which allowed the assessment of the specificity of ONE-seq platform. Other samples that contained 1%, 5%, and 10% of NAD-relative to m 7 G-RNA spike-ins were used to determine the capture sensitivity. We subjected 100 g total RNA from mouse livers mixed with 1 ng spike-in RNA to ONE-seq experiment, followed by polyA-selected RNA sequencing. In the sample mixed with 100% m 7 G-RNA spike-in, we found no enrichment ( Figure  3I). As demonstrated by sequence read count, the amount of spike-in RNA exceeded 99% of endogenous mRNA transcripts from mouse livers (Supplementary Figure 4A), representing an abundant transcript in the transcriptome. This evidence highlights the specificity of ONE-seq, which can eliminate potential contamination from m 7 G-capped transcripts that are highly expressed. In the sample that contained 1% of NAD-capped forms, the enrichment was low and variable ( Figure 3I), suggesting that ONE-seq might not be sensitive enough to capture low-degree NAD modification for particular transcripts. In contrast, when NADcapped forms accounted for 5% or 10% of the spike-in transcript, the level increased up to 3.4-fold and 4-fold, respectively ( Figure 3I), reflective of significant enrichment. This experiment provided an estimate of the stoichiometry of NAD versus m 7 G in the endogenous transcripts by leveraging spike-in RNAs with ascending ratios of NAD-capped forms. As a result, we proceeded to set 2-fold enrichment, roughly reflecting 3% of NAD-capped form for a particular transcript, between ONE-seq and input as the cutoff for the identification of NAD-RNAs.
Fourth, we determined the noise-cancelling effect of NudC by comparative analysis of RNA-seq datasets. Total RNAs, after HEEB reaction, were captured by streptavidin beads, followed by either NudC-catalyzed elution (NudC+) or mock treatment (NudC-) ( Supplementary Figure S3). polyA-selected RNA sequencing was performed. To pinpoint NAD-capped RNA, we set 2-fold enrichment of read counts as the cutoff (Supplementary Figure S3). From mouse liver tissues, 1,952 NAD-RNAs were identified in the absence of NudC; however, the usage of NudC expelled 489 genes from the profile ( Figure 3J and Supplementary Table S1). As a result, NudC treatment may contribute by roughly reducing 25% of the false positive hits presumably derived from m 7 G-capped RNAs. (G) qRT-PCR analysis shows that NAD-RNA (106 nt), but not ppp-RNA (106 nt), can be enriched by ONE-seq. (H) qRT-PCR analysis shows that streptavidin beads bound HEEB-reacted m 7 G-capped RNA (106 nt) cannot be eluted by NudC (Two-tailed Student's t test: ***P < 0.001; **P < 0.01; n.s., not significant). (I) RNA-seq experiment of spike-in RNAs determines the sensitivity of ONE-seq. Top panel: schematic workflow of total RNAs and polyadenylated spike-in RNAs that had different ratios of NAD-RNA (500 nt). Two spike-in RNAs with identical sequence (500 nt) but have either NAD or m 7 G-cap, followed by polyA tails are used; bottom panel: fold change of normalized read counts from spike-in RNA between enrichment and input samples in different ratios of NAD-RNA. Total RNAs were from liver tissues of 18-month mice. The nominal ratios of NAD-RNA were highlighted in blue. (J) Epitranscriptome assessment of NudC to minimize the noise of m 7 G-RNAs. Two-fold enrichment of read counts was used as the cutoff. Standard ONE-seq identified 1,638 NAD-RNAs, while 489 false-positive NAD-RNAs were found without the use of NudC, presumably derived from m 7 G-capped RNAs. Total RNAs were from liver tissues of 12-month mice.

Epitranscriptome-wide analysis of NAD-RNAs by ONE-seq
We tested the utility of ONE-seq. We profiled NAD-RNAs from mouse liver tissues of young (2-month) and aged (18-month) cohorts. After quality control, we obtained in average ∼12.3 million high-quality and uniquely mapped sequencing read pairs from each library (Supplementary Figure S4B). Assessment of datasets corroborated that sequencing saturation has been reached (Supplementary Figure S4C). Principal component analysis (PCA) illustrated that the biological replicates of each sample type clustered together, reflecting high reproducibility of the experiments (Supplementary Figure S4D). We proceeded to set 2-fold enrichment of read counts as the cutoff, which led us to identify 2017 and 1820 NAD-RNAs from young and aged animals, respectively (Supplementary Table S2). Notably, similar distributions of read counts along gene bodies were shown between mRNA transcriptome (input) and NAD epitranscriptome (ONE-seq) (Supplementary Figure S4E), suggesting comparable integrity of transcripts between m 7 G and NAD-capping.

Validation of NAD-RNAs by boronic acid
We validated newly identified NAD-RNAs by an ADPRCindependent method. To do this, we applied boronic acid beads. Our rationale lies at the fact that boronyl groups of boronic acid can form relatively stable complexes with 1,2-cis diols, occurring naturally at the 3 -end of RNA, as well as in the 7-methylguanosine of m 7 G-cap and the nicotinamide riboside of NAD-cap (37)(38)(39). Inspired by previous assays (37,38), we devised a qRT-PCR strategy to specifically capture NAD-RNA ( Figure 4A). First, yD-cpS is used to remove 7-methylguanosine at the 5 -end of m 7 G-RNA, while nicotinamide riboside of NAD-cap remains intact. Second, an RNA adaptor that contains a 2 ,3dideoxycytidine (ddC) at the terminal is ligated to the 3end of RNA, such that the vicinal-diol moiety of the ribose at the 3 -end of RNA no longer exists. At this step, affinity binding can only occur between the boronyl group from boronic acid and 1,2-cis diols from the nicotinamide riboside of NAD-cap. Importantly, we perform reverse transcription using an adaptor-specific primer to harvest RNAs with 3 end blocker, thereby abolishing false signals from boronate binding of RNAs that are not properly ligated with adaptor. Consequently, only RNA transcripts that contain NAD-cap can be selectively enriched by boronate affinity.
As shown in Supplementary Figure S2A, yDcpS treatment specifically de-capped m 7 G, but not NAD, from RNA oligos (38 nt). In the presence of RNA ligase, 3 adaptor was efficiently ligated to the spike-in RNAs (38 nt), yielding an upper band in the gel ( Figure 4B). Upon these treatments, we showed the evidence that NAD-capped, but not m 7 G-capped RNA oligos, were retained by boronic acid beads ( Figure 4C). We then extended this strategy to the total RNA extracts from mouse livers mixed with three long spike-in RNAs (ppp-RNA and m 7 Gppp-RNA as negative controls, NAD-RNA as positive control) that were designed with different sequences. This experiment clearly demonstrated that NAD-RNA (500 nt) spike-ins, but not ppp-RNA (500 nt) and m 7 Gppp-RNA (500 nt), were selectively and significantly enriched by boronic acid beads (Figure 4D). We proceeded to successively validate the NADcapping for four endogenous transcripts, Akr1c13, Atp5k, Prxl2b, and Med17, as identified by ONE-seq ( Figure 4D). In contrast, two genes that were not identified by ONE-seq, Serpinale and Fgb, showed no enrichment by boronic acid beads ( Figure 4D). Combined, this ADPRC-independent validation corroborated newly identified NAD-RNAs, supporting the notion that ONE-seq, an ADPRC-dependent method, enables robust and efficient identification of NAD-RNAs.

Characterization of NAD-RNAs from mouse livers
We characterized newly identified NAD-RNAs from mouse livers by ONE-seq. We extracted top 100 most abundant transcripts from input dataset and examined their relative enrichment. This analysis found no evidence of positive correlation, indicating that NAD capping might not be simply incidental events that are proportionally correlated with transcript abundance ( Figure 5A). In mouse livers, NAD capping mostly occurred on protein-encoding genes, but also extended to non-coding RNAs and pseudogenes ( Figure 5B and C). NAD-RNAs were shown to be derived from genes localized on autosomes and X chromosomes, but not from the Y chromosome and the mitochondrion genome ( Figure 5D). By dividing NAD-RNAs into 10 deciles based on enrichment, we observed that shorter genes tended to have increased modification of NAD ( Figure 5E). Analysis of the splicing events revealed a higher ratio of intron retention in NAD-capped than those of non-NAD-capped mRNA transcripts ( Figure 5F). Genome browser views illustrated the presence of intron read counts in mRNAs capped by NAD, whereas those introns were normally spliced in non-NAD capped forms ( Figure 5G).
To gain functional insight, we performed pathways analysis on RNAs found to contain NAD-cap. We started by analyzing dataset from young mouse livers. Analysis of biological processes and molecular pathways revealed that NAD-RNAs were mainly involved in DNA replication, transcription, cell cycle, metabolism, immune system, response to stimulus, and ribosome biogenesis ( Figure 6A). In addition, genes linked to subcellular localization and function, such as those in nucleus and mitochondrion, were highly enriched.
In nucleus, genes encoded basic machineries in DNA replication (e.g. DNA primase and DNA polymerase) and transcription (e.g. RNA polymerase II) were found to contain NAD-cap. Also, transcripts with 5 -NAD modification were enriched in cell cycle (e.g. centromere protein and CDK) and DNA damage responses (e.g. ATM and TP53). Not just housekeeping genes, transcriptional factors (e.g. BMP5 and Interferon), epigenetic pathways (e.g. Mediator Complex, DPY-30 and HMGN2), and histone-related genes were also subject to NAD modification ( Figure 6A and Supplementary Table S3). In mitochondrion, we noted that its basic transcription (e.g. mitochondrial transcrip- In the presence of yDcpS, 7-methylguanosine of the m 7 G-cap was removed. Prior to 3 -end ligation, RNAs were treated with CIAP to remove 5 -terminal phosphate. Adapters with 2 ,3 -dideoxycytidine (ddC) were ligated for blockage of the cis-diols moiety at the 3 -end of RNA. At this step, affinity binding can only occur between the boronyl group from boronic acid and 1,2-cis diols from the nicotinamide riboside of NAD-cap. Adaptor specific reverse transcription, followed by gene-specific qPCR, was performed to assess the NAD modification. In the right panel, 1,2-cis diols were highlighted in red dash rectangles. (B) Ligation of 3 adaptor resulted in the appearance of upper bands. (C) NAD-RNA, but not m 7 Gppp-RNA RNA oligos, were retained by boronic acid beads, after yDcpS treatment and 3 -adaptor ligation. (D) Assessment of gene-specific NAD-capping by qRT-PCR. Based on boronic acid beads, Akr1c13, Atp5k, Prxl2b, and Med17 identified by ONE-seq as well as Serpinale and Fgb not identified by ONE-seq were examined. (Two-tailed Student's t test: ***P < 0.001, *P < 0.01; n.s., not significant). tion termination factors) and translation (e.g. mitochondrial ribosomal proteins) were highlighted by NAD-RNAs. Nuclear-encoded mitochondria genes, known to play essential roles in electron transport (e.g. mitochondrial complex I and cytochrome P450) and metabolism (e.g. SDHC and LDHB), had NAD modification ( Figure 6A and Supplementary Table S3). In contrast, no transcripts encoded by the mitochondrial genome were found to be linked with NAD.

Age alters NAD-capped RNAs
In line with an age-modulated decrease of NAD (Supplementary Figure S5), both the number of RNAs that contained NAD-capped form and the global extent of capping became decreased in aged compared to young animals (Figure 6B). Gene-specific analysis conformed with this trend of age-modulated decrease in NAD capping; however, select NAD-RNAs had higher enrichment in aged than young mouse livers ( Figure 6C). Interestingly, we noted different and 1820 NAD-RNAs from young (2-month) and aged (18-month) mouse livers, respectively. Total RNAs were from mouse livers of indicated age. (C) NAD capping on different RNA types, with most occurrence on protein-encoding genes (blue), but also on non-coding RNAs (orange) and pseudogenes (yellow). (D) Chromosomal distribution shows that NAD-RNAs are derived from genes localized on autosomes and X chromosomes, but not from the Y chromosome and the mitochondrion genome. (E) From 10 deciles based on enrichment, genes with short length tend to have increased modification of NAD. (F) NAD-RNA tends to have higher ratio of intron retention than non-NAD capped forms. (G) Genome browser views illustrate the presence of intron read counts in select NAD-capped genes.  functional categories in young and aged animals. Whereas young animals displayed enrichment of pathways involved in nucleotide excision repair (e.g. excision repair protein) and RNA splicing (e.g. splicing factor and alternative splicing regulator), aged animals tended to enrich metabolic signatures, such as ribonucleotide metabolism (e.g. UPP2) and aromatic-amino acids metabolism (e.g. TDO2), as well as oxidation-related pathways (e.g. GSTs and SODs) (Figure 6D and Supplementary Table S3). Together, our study reveals new features of age-modulated NAD-RNAs from adult mouse livers.

ONE-seq enables gene-specific analysis of NAD-RNAs
Above evidence supported the notion that ONE-seq platform permits epi-transcriptome-wide profiling directly from total RNA, prompting us to extend its application for genespecific assessment by qRT-PCR. To do this, we included non-capped ppp-RNA (106 nt) as a baseline negative control, and NAD-RNA (106 nt) as a positive control. Total RNA and internal spike-in controls were subjected to ONEseq experiment, followed by qRT-PCR. The relative abundance was calculated between NudC + and input samples ( Figure 7A). Our data demonstrated that ONE-seq enabled comparative and quantitative assessment of NAD-capping events on specific genes, such as the cytochrome P450 family Cyp2c70 and Cyp3a11 involved in electron transport, Akr1c13 and Prxl2b of the metabolism-related genes, as well as Med17 and Ufc1 genes of gene regulatory pathways ( Figure 7B).

DISCUSSION
NCIN, the nucleoside-containing metabolite such as NAD, 5 -desphospho coenzyme A (dpCoA), flavin adenine dinucleotide (FAD), uridine diphosphate glucose (UDP-Glc) and uridine diphosphate N-acetylglucosamine (UDP-GlcNAc), can be incorporated at the 5 -end of RNA during transcription initiation in both prokaryotes and eukaryotes (5,34). CapZyme-seq has been applied to detect the global landscape of NCIN-capped RNAs in vivo but this method does not pinpoint the nature NCIN cap (34). As a result, CapZyme-seq cannot be used to profile NAD-RNAs from total RNA extract.
NAD is the hub metabolite and redox regent for cells, involving in a wide range of biological processes (9). The attachment of NAD to RNA inherently connects essential metabolic regulation with gene expression, defining a critical layer of epi-transcriptomic regulation. However, investigating biological insight of NAD-RNAs has been hindered by the analytical methods available. The currently reported NAD-RNA identification methods involve the use of multiple reactions, and each reaction requires additional steps of RNA clean-up and precipitation. In addition, the canonical 5 -end cap structure of RNA (m 7 G) has been found to contaminate the NAD-RNA profile. As such, SPAAC-NADseq (16), the most recent method developed in Arabidopsis, introduces antibody-based pre-treatment to deplete m 7 G-RNA from purified mRNA. Consequently, these methods demand laborious procedures and high RNA input, which cannot be readily applied for gene-specific analysis. In the current study, we design a one-step chemo-enzymatic reaction by HEEB that directly conjugates NAD-RNA with biotin affinity tag. We apply NudC-based post-treatment to specifically harvest NAD-capped RNAs, allowing the assay to be simply performed from total RNA. Compared to previous methods, ONE-seq requires significantly less amount of RNA input. More importantly, we include different types of spike-in RNAs to carefully demonstrate that ONE-seq can capture NAD-RNAs with specificity and sensitivity, especially for those that have equal to or higher than 3% of NAD-capped form for a particular transcript (see Figure  3I). However, one limitation to note is that ONE-seq is not sensitive enough to detect low-degree NAD modification, such as those equal to or lower than 1% of the transcript (see Figure 3I). For future experiments, synthetic m 7 Gand NAD-RNA spike-ins can be included as quality controls to determine specificity and sensitivity, respectively. Inspired by previous work (38,39), we devise an ADPRCindependent assay by boronate affinity to validate NAD-